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Abstract 

The invariance theory in continuum mechanics is applied to analyze 
Reynolds stresses in high Reynolds number turbulent flows. The anal- 
ysis leads to a turbulent constitutive relation that relates the Reynolds 
stresses to the mean velocity gradients in a more general form in which 
the classical isotropic eddy viscosity model is just the linear approximation 
of the general form. On the basis of realizability analysis, a set of model 
coefficients are obtained which are functions of the time scale ratios of the 
turbulence to the mean strain rate and the mean rotation rate. These 
coefficients will ensure the positivity of each component of the turbulent 
kinetic energy — realizability that most existing turbulence models fail to 
satisfy. Separated flows over backward-facing step configurations are taken 
as applications. The calculations are performed with a conservative finite- 
volume method. Grid-independent and numerical diffusion-free solutions 
are obtained by using differencing schemes of second-order accuracy on suf- 
ficiently fine grids. The calculated results are compared in detail with the 
experimental data for both mean and turbulent quantities. The compar- 
ison shows that the present proposal significantly improves the predictive 
capability of K-e based two equation models. In addition, the proposed 
model is able to simulate rotational homogeneous shear flows with large 
rotation rates which all conventional eddy viscosity models fail to simulate. 
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1. Introduction 


Numerical simulation of turbulence is a bottleneck in the development of compu- 
tational fluid dynamics (CFD). The approach of direct numerical simulation (DNS) 
without any turbulence models is restricted to simple and low Reynolds number 
flows within the capabilities of current computers. A compromise to DNS is the 
large-eddy simulation (LES) approach in which the large scales of turbulence are di- 
rectly computed and the small scales are modeled. LES, though applicable for high 
Reynolds number turbulence, is usually very expensive and has serious problems 
with boundary conditions. Some of these difficulties also exist in DNS. Therefore, 
most practical calculations at the present time are based on averaged Navier-Stokes 
equations with the aid of turbulence modeling. 

In turbulence modeling, unknown turbulent correlations are expressed in terms 
of determinable flow quantities with the aid of empirical information. According to 
the way the Reynolds stresses (the second-order moment correlations) are treated, 
turbulence models may be divided into two groups: the Reynolds stress algebraic 
equation models and the Reynolds stress transport equation models. The former 
group includes the zero-, one- and two-equation models (Rodi, 1980) in which the 
Reynolds stresses are algebraically related to the mean flow field. The eddy vis- 
cosity K-c two-equation model (Launder and Spalding, 1974) in this group is one 
of the most popular turbulence models used today in practical flow calculations. 
In the latter group, often called second-order closure models, the Reynolds stresses 
are determined by their own dynamical transport equations. Second-order closures 
are attractive because they can simulate the transport of the individual Reynolds 
stresses; however, it is difficult to consistently model all the higher-order turbulent 
correlations appearing in these second moment dynamical equations. Inappropriate 
modeling of higher order correlations (often due to lack of information about their 
underlying mechanism) could result in a serious inaccuracy and unphysical results. 

In the standard K-c model, all the model coefficients are constant and are de- 
termined from a set of experiments for simple flows under equilibrium or isotropic 
turbulence conditions (Rodi, 1980). Numerical experience over the last two decades 
has shown that this set of constants has a broad applicability, but this by no means 
signifies that they are universal. Rodi (1980) found that the K-e model’s ability to 
predict weak shear flows can be significantly improved by using as a function of 
the average ratio of P/t (P is the production of the turbulent kinetic energy) instead 
of a constant. Leschziner and Rodi (1981) proposed a function for C „ which takes 
into account the effect of streamline curvature and obtained improved results in the 
calculation of annular and twin parallel jets. Recently, Yakhot and co-workers have 
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developed a version of the K-e model using renormalization group (RNG) methods. 
Their model is of the same form as the standard K-e model, but all the model coef- 
ficients take different constant values. In the latest version of the RNG based K-e 
model (Speziale and Thangam, 1992), the coefficient Ci, related to the production 
of dissipation term, is set to be a function of 77 , where 7 ) is the time scale ratio 
of the turbulent to mean strain rate. In applying this model to a separated flow 
over a backward-facing step, experiment ally studied by Kim et a J. (1978), Speziale 
and Thangam obtained a good prediction of the reattachment length which is an 
important parameter often used to assess the overall accuracy of calculations. 

The standard K-e model (including the RNG based one), like many others in the 
algebraic equation model group, uses Boussinesq’s isotropic eddy- viscosity concept 
which assumes that the Reynolds stresses are proportional to the mean velocity 
gradients. The concept usually does well for the shear stresses in two-dimensional 
mean flows of the boundary-layer type, but not well for the normal stresses due to 
the erroneous isotropic nature of the concept. This suggests that linear dependence 
on the mean velocity gradients is insufficient and that a more general relation is 
needed for more complex flows. In fact, by eliminating the convection and diffusion 
terms in the modeled transport equations for the Reynolds stresses, Rodi (1980) 
developed an algebraic stress model (ASM) in which the Reynolds stresses are cal- 
culated by algebraic expressions. Owing to its anisotropic nature, the model does 
perform better than the isotropic K-e model for certain flows; a well known example 
is fully-developed flow in non-circular ducts where ASM is capable of generating 
turbulence-driven secondary motions while the isotropic eddy viscosity K-e model is 
not. However, ASM does not appear in a tensorial invariant form, which may limit 
its generality. In addition, inappropriate modeling of higher order correlations, such 
as pressure-strain correlations, will also cause deficiencies of the second-order clo- 
sure based ASM. Moreover, special care needs to be taken to prevent the turbulent 
normal stresses from becoming negative (Huang and Leschziner, 1985), and the nu- 
merical implementation of ASM may even be more complicated than that of its 
parent second-order closure model, especially in general three dimensional flows. 
Recently, this numerical difficulties of ASM was first nicely resolved by Taulbee 
(1992). 

There are other approaches to developing Reynolds stress algebraic equation mod- 
els. For example, Yoshizawa (1984) derived a relation for the turbulent stresses using 
a two-scale direct interaction approximation. It contains both linear and quadratic 
terms of the mean velocity gradients. A similar relation wets also derived recently by 
Rubinstein and Barton (1990) using Yokhot and Orszag’s RNG method. An inter- 
esting point in these two methods is that the values of the model coefficients can all 


3 



be determined analytically. Speziale (1987) proposed a different expression, based 
on the principle of material frame-indifference, which contains the Oldroyd deriva- 
tive of the mean strain rates. However, the principle of material frame-indifference is 
only valid in the limit of two-dimensional incompressible turbulence, hence it is not 
an appropriate constraint for general turbulent flows. In addition, these non-linear 
models are not fully realizable and have not been extensively tested. 

The purpose of the present study is to develop a general and realizable Reynolds 
stress algebraic equation model with the method of rational mechanics. As usual, 
we assume that the Reynolds stresses depend on the mean velocity gradients, the 
turbulent velocity and length scales, then a constitutive relation for the Reynolds 
stresses is derived by using the invariance theory. The final form is truncated up 
to tensorial quadratic terms of the mean velocity gradients. Using the re aliz ability 
conditions, the coefficients in the obtained relation are found to be at least functions 
of the time scale ratio of the turbulence to the mean strain rate. In general, they 
are also functions of the time scale ratio of the turbulence to the mean rotation rate. 

The model validation is made on the basis of applications to the rotational ho- 
mogeneous shear flows simulated by Bardina et a 1. (1983) and the two backward- 
facing step flows experimentally studied by Driver and Seegmiller (1985) and Kim 
et ai. (1978). The latter type of flows has served as a benchmark in validating 
turbulence models for complex flows. Calculations are carried out with a conser- 
vative finite-volume method, and a second-order accurate and bounded differencing 
scheme, together with sufficiently fine grids, is used to ensure that the solution is 
both grid-independent and free from numerical diffusion. The calculated results 
are compared in detail with experimental data as well as with those obtained using 
standard K-e model. 


2. Modeling of Reynolds Stresses 

Incompressible turbulent flows are governed by the following Reynolds averaged 
continuity and Navier-Stokes equations: 

Ui ti = 0 ( 1 ) 

Ui,t + (UjUi - vUij -I- uiuj)j = (2) 
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where U% are the mean velocity components (t = 1,2,3), p is the mean pressure, v 
and p are the fluid kinematic viscosity and density, Ui,t and Uij are derivatives of Ui 
with respect to time t and co-ordinate x,j, respectively, uiuj is the turbulent stress 
tensor which must be modeled. 

The oldest and simplest proposal for modeling the turbulent stress is Boussi- 
nesq’s isotropic eddy- viscosity concept that assumes an analogy between the viscous 
stresses in laminar flows and the turbulent stresses in turbulent flows. The general 
form of this concept is 

iuuj = — ^ v tSij (3) 

where »/ t is called the eddy-viscosity and Sij is the mean strain rate defined by 

% = \(Vij + U,,) (4) 

Equation ( 3) constitutes a common basis for most turbulence models that are 
extensively used today. 

2.1 Constitutive relation 

Does a general constitutive relation exist for turbulent correlations? Lumley 
(1970) discussed this problem and found that such relations exist only under the 
situation in which the length and time scales of turbulence are much smaller than 
those in the mean flow field so that the effect of initial and boundary conditions 
on the turbulence is not significant far from the wall. In other situations such as 
rapidly developing mean flows or in the vicinity of walls, it is questionable whether 
there exists such a constitutive relation for any turbulent correlation; however, from 
practical point of view, we can formally derive a “constitutive” relation for any 
turbulent correlation to solve the closure problem. The validity of such a formally 
derived relation needs, of course, to be verified with the aid of experiments. 

A turbulent constitutive relation, if it exists, is always of functional form. From 
a modeling point of view and for convenience of application, we neglect the time 
memory effects in the relation and consider the relationship at the present time as the 
first order approximation in the time expansion of the functional form. Therefore, we 
assume that the turbulent stress U{Uj is a function of the mean deformation tensor 
Uij , the velocity and length scales of turbulence characterized by the turbulent 
kinetic energy K and its dissipation rate e, i.e., 

u^Uj = Fij(Uij, K, e) (5) 
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Note that the molecular viscosity v is not included because we restrict our attention 
here only to high Reynolds number turbulent flows. 

The arguments of equation (5) contain ten quantities bearing two dimensions. 
According to the ic theorem of dimensional analysis, they may be grouped into eight 
independent non-dimensional quantities: 

= jVij ( 6 ) 

By normalizing the turbulent stress as 


equation ( 5) can be written as 


ViVj 


UjUj 

2K 


( 7 ) 


m = Wa) ( 8 ) 

The form of the tensor valued isotropic function F,j can be determined by using 
invariance theory (Lumley, 1978). The basic principle is that an invariant can only 
be a function of other invariants. In determining a set of independent invariants, we 
have shown in Appendices A and B (also see Shih, 1992) that only 18 independent 
tensors can be formed with the tensor Vij and its transpose Vjj according to the 
generalized Cayley- Hamilton relations. Following Lumley (1978), let A, and Bj be 
two non-dimensional arbitrary vectors, we may form the following 18 invariants: 


VijAiBj t VjjAiBj, 

V iik V hij AiB iy Vj'kVk'iAiBj, Vi'kVj'kAiBj, V hti V k jAiB h 
V itk V/ tk AB jf V? k Vj t kAiBj, Vk'iVfjAiBj, VfrVLjAiB,, 
V i *V,, k V,' J A.B h V&VfrMBi, V&VhJUB,, 
Vvy&V&AiB,, 

Vt.kV,j,V,'„V j * m A i Bj, V?,v;,V„jV mj A,B, 


( 9 ) 


In addition, we have other invariants: 


SijAiBj, AiAi , BiB iy /, II, ///, 


where I, II, III are the three invariants of the tensor Vij: 


( 10 ) 
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( 11 ) 


/ = V iti 

u = \iViMj - VijV^i) 

III = ^(ViMjVk'k - 3 ViMtVkj + 2 VijVijMs) 

and • • • represents the invariants of other 17 tensors, for example, VijVij. 

Note that in the above list of invariants we do not include Vi t iVi t kVi,jAiBj and 
higher order terms of this type because they are not independent quantities according 
to the Cayley-Hamilton theorem: 


V„V,. k V kJ - I ■ V,j.V kJ + II ■ K., - Ills,, = 0 (12) 

Any other possible terms, for example, V kl iV k jV mi iV m jAiBj, are also not indepen- 
dent, hence they will not be included. 

However, the invariant list can be extended by including any combination of 
invariants in ( 9) and ( 10), for example, 


VijAiBj-MI, II, III , V^, .••) 


(13) 


Vi'MjAiBj ■ f 2 (I, II, III , V^, •••) 

where f\ and fi are scalar functions. Of course, these types of invariants axe not in- 
dependent, but they are useful in explaining why the coefficients in the final relation 
( 16) are, in general, functions of various invariants of the tensors in question. 

As a result, the invariant v^VjA,Bj may be written as a function of the above 
invariants listed in ( 9), ( 10) and ( 13): 


ViVjAiBj = /( VijAiBj , VjjAiBj, Vi t kVkjAiBj, Vj^Vk^AiBj, 
Vi+VjkAiBj, VkjVkjAiBj, Vi, k V/ ik AiBj, 

V^AiBj, V k>i \ tfjAiBj, Vt-VujAiB;, 

Vi'kVi'kVfjAiBj, V? iV^AiBi, VltfjAiB), 
Vk.iVk'iVfjAiBj, Vki Vfj V*j Ai Bj , V^V^V^jAi Bj , 
Vi,kVi,kV* m Vj m AiBj, VlM tl V m jVnjAiBj, 

SijAiBj, AiAi, B,B, , I, II, III , VijVij, •••, 
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V,.,A,B,f u ...) 


(14) 


Because UJtJJ AjBj is bilinear in arbitrary tensors A{ and Bj y we must require that 
the right hand side of equation ( 14) be also bilinear in Ai and Bj. Therefore, 
equation ( 14) can be reduced to 


ViVjAiBj — a j 8 y j A y Bj -J- n j A( Bj -|- Vj t % A+ Bj 

+ <*5 Vj'kVk'iAiBj + atVijtVj'kAiBj 

+a r V k , i V kJ A i Bj + + a^V^ABj 

+a I „V Mi V!jAB i + « nVbVuAiBi + ^V iJ .V, J .V^A i B j 
+^V'kV^A,B, + a^V’jA.B, + *„V i .<V i jV? J A.Bj 
+a l ,V ii ,V’jV; t A i B i + aMtftfjAiBi 

+^.V i , k V l ,,.V t ] m Vf m A i B i + « t ,V i ’ i V^V„ J V mJ A i B i (15) 

Noting that Ai and Bj are the arbitrary vectors, we obtain 


V i V j — a l &ij + diViJ + ®3 Vj,i + &<Vi,kVkj + 0-f,Vj ik V k j + 0-eVi >k Vj t l k 

+a 7 V k ,iV kii + a*V iM Vl k + a 9 V? k V iik + a 10 V k>i V^ + 
+a 12 V itk V l , k V l i j + a 13 V* k V? k + a lA V^. + 

+ai 6 V kii V* l V* l + &irVi tk V* k V*j + 

+ a uV k ,i V kt i V m j V m j 


(16) 


where the coefficients aj — are, due to ( 13), functions of the invariants /, //, 
III, VijVij • • •, i.e., 


<*< = fi(I, II, III, VijVij , • • •), * = 1 , 2 , — , 10 

Equation ( 16) is a general relationship between two second rank tensors. 


(17) 
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The normalized turbulent stresses Viv] have two properties 


and 


W>j = VjVi 


( 18 ) 


ViVi = 1 (19) 

Using the properties of ( 18) and ( 19) in equation ( 16), we obtain the following 
relations: 


“3 — 




CL\2 = <Ii6 = ®1S — a 17 — a 18 — °19 — 9 

oi = —[1 — 2aj I — 2a 4 D — (ae + o-i )D — 2(a# + ujo )D — (uis + oujft] 

O 


( 20 ) 


where 

D = V U V M , D = V.jV.j , 

After introducing equations ( 20) into ( 16) 
form, we obtain 


b = VijV?j , D = V?jV?j (21) 
and converting to the dimensional 


UiU j 


\K6„ + 2 a,~(U u + V u - jDiAi) 
+2 + U iM Vu,i - |n«iy) 


K 3 1 

+2a,— (Ui.kUjj, - -US.,) 


1 - 


K 3 1 

+2a 7 —(V k , i U>j - 
+2c a ^-(U i . i Ul + DJjOm - 
+2a,„ ^(V ki Ul 4 + VkjVli ~ 
+2«„^(l 


where 


(22) 


+2 - ^ft Sa) 

n = UijUj'i , n - u^j , 6 = ^ 5 , ft = u^ul (23) 
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Equation ( 22) is the most general relationship between and U%j under the 
assumption ( 5). Interestingly enough, the first five terms at the right-hand side 
of equation ( 22) are of the same form as those derived through both the two- 
scale direct interaction formalism (Yoshizawa, 1984) and the renormalization group 
method (Rubinstein and Barton, 1990). The fact that the three different theoretical 
analyses lead to a similar result indicates the rationality of equation ( 22). 

In practice, however, a quadratic tensorial form of ( 22) may be sufficient, espe- 
cially when \\Uij\\K/e is less than unity (which is true if the turbulent time scale is 
smaller than the time scale of mean flow). Therefore, from now on, we only consider 
the quadratic form of equation ( 22). 

2.2 Realizability 


Realizability (Schumann, 1977, Lumley, 1978), defined as the requirement of 
the non-negativity of turbulent normal stresses and Schwarz’ inequality between 
any fluctuating quantities, is a basic physical and mathematical principle that the 
solution of any turbulence model equation should obey. It also represents a minimal 
requirement to prevent a turbulence model from producing unphysical results. In 
the following, this principle will be applied to the constitutive relation ( 22) to derive 
constraints on its coefficients. 

Consider a deformation rate tensor of the form 

/ U XA 0 0 \ 

0 U 2t2 0 (24) 

0 0 0 / 

The continuity equation (1) gives 


U 2t2 = —U 2f i 

and from equation ( 22), the normal stress UiU\ can be written as 


(25) 


UiUj 

~1K 


= i + + i(2a, + a. + a,) 


(26) 


Since the time scale ratio of the turbulent to the mean strain rate is defined by 

KS 


V = 


(27) 


where 


S = (2 SySy) 1 '* 


(28) 
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equation ( 26) can further be written as 


~Tir = 3 + air l + I2^ 2a4 + a# + a 7 )^* ( 2Q ) 

Physically, we know that U\U{ will decrease by a vortex stretching with an increase 
in U\ t i, but uiUi cannot be driven to negative values. Therefore, we must require 
that 


«itti 

2K 

uTuT 

2K 


\ 2K ),n 


> 0, 

if 0 < 17 < 00 

(30) 

- 0, 

if 77 —> 00 

(31) 

- 0, 

if 77 — ♦ 00 

(32) 


These are called the realizability conditions. They can be satisfied in various ways 
of which the simplest way is perhaps the following: 


2a 3 — 

2a 4 — 

2a e = 

2a 7 = 


2/3 

A i+i] 
Ct\ 

f(v) 

c r2 

f(v) 

CrZ 

m 


(33) 

(34) 

(35) 

(36) 


where f(rj) is in general a polynomial of rf of order higher than 2. We take its 
simplest form as 


f(v) = A 2 + v s (37) 

Ai, A 2 ,C r i,C T 2 and C T 3 are adjustable constants, but they must satisfy 


A-i > 0 , A 2 > 0 , 

2C T 1 + C t2 + CVs > o 


(38) 


Similar analysis on U 2 U 2 and u 3U3 also leads to equations ( 33) to ( 38). It should 
be mentioned that equations ( 33) to ( 38) also hold for a three-dimensional pure 
strain rate tensor 
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( 39 ) 


l U IA 0 0 \ 

0 U 3>3 0 

V 0 o u 3 , s J 

and that any deformation rate tensor can be written in the form of ( 39) in the 
principal axes of deformation rate tensor. 

It can be seen from the above analysis that realizability cannot be fully satisfied 
if the model coefficients are taken as constant, such as those in the standard K-e 
model and in the recent anisotropic models of Speziale (1987), Yoshizawa (1984) 
and Rubinstein and Barton (1990). In fact, these models satisfy realizability only 
in the weak sense, that is, they only ensure the positivity of the sum of the normal 
Reynolds stresses. 

Further constraints on the model coefficients can be obtained by considering the 
deformation rate tensor 


0 U X j 0 
0 0 0 
0 0 0 

which corresponds to a fully- developed channel flow. In this case, we have 


(40) 


where 


2 t] 2 K 

U lUl - -K + ^^^(^ " Cr 3 ) 


U3U3 

UiU 2 


l K £k 

3 3(A 2 + t/3) 

2r}K 


(CtJ + C T 3) 


3(j4i -+■ 77 ) 


V = 


KS 


S = (2 = |0i 4 | 


Experiments indicate that 


(41) 


(42) 


3 

2 „ 

u 2 u 2 < -K 


which requires, from equation ( 41), that 


(43) 
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C r 2 > 2C r 3 


(44) 


2.3 Rotation effect 


The parameter 77 represents the effect of the mean strain rate, and the effect of 
the mean rotation rate can be represented by 


$ = n = fiv = (v u - f/,,)/2 + 4 £ „„u, m (45) 

where u> m represents the rotation of the frame. 

In the present study, we find that it is sufficient to simply include the parameter 

£ only in the coefficient a2, i.e., 

2a, = -- 2/3 ' - («) 

Ai+tj + a£ 

while keeping the other coefficients the same. The dependance of the coefficients on 
17 and £ can be easily justified by equation ( 17). 


2a? — — 


2.4 Realizable algebraic equation model 

By introducing equation ( 33)-( 37) into equation ( 22), we obtain 


muj = -K6ij — Vt{Uij + Uj,i) 

+ (47 . 

A 2 + TJ e 2 3 

+t - ^n<ij) 

M + T e 6 

Two quantities, the turbulent kinetic energy K and its dissipation rate e, remain to 
be determined in equation ( 47). To this end, the two transport equations in the 
standard K-e model are used which read: 


K, t + [UjK ~(v+ —)Kj]j = P-e 

&K 


tjt + Wit -(v+ j = c *r p - ° 2 V 
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where 


( 50 ) 


£1 2/3 

t ’ ' A,+.| + a( 

P = -Oju jU,j (51) 

The coefficients C \ , C ? , &K and c t assume their standard values: 

Ci = 1.44, C 2 = 1.92, <t k = 1, <r t = 1.3 (52) 

and the additional coefficients assume: 

Cri = -4, Cr2 = 13, C r3 = -2, Ar = 1.25, a = 0.9, A 2 = 1000. (53) 

These values have been found to work well for both test cases considered in this 
work. 


3. Rotating Homogeneous Shear Flow 

The present model is able to mimic the effect of the mean rotation rate on the 
turbulence. A test case is the rotating homogeneous shear flow which was studied 
by Bardina et ai. (1983) using the LES method. Figure 1 is the configuration of 
the flow being tested. Figures 2(a - c) show the evolution of the turbulent kinetic 
energy K/Ko with the nondimensional time, St, at the rotation rates of Cl/S = 
0, 0.5, —0.5 respectively, where Ko is the initial turbulent kinetic energy, 5 is the 
mean strain rate and 17 is the rotation rate of the reference frame. The calculations 
were performed with a fourth order Runge-Kutta scheme. The initial condition 
corresponding to the isotropic turbulence used in LES with cq/(S K o) = 0.296 was 
adopted for all the three cases. The results from both the present model and the 
standard K-e model (hereafter referred to as s-K-e) are compared with LES results 
in figures 2(a - c). It can be seen that at Q/S = 0 the present model cannot 
predict the initial nonequilibrium development of turbulent kinetic energy very well 
starting from an isotropic turbulence. However, it does catch up with the later 
“equilibrium” development and performs much better than the s-K-e model which 
highly overpredicts the data. Figures 2(b) and 2(c) show the ability of the present 
model to simulate the effect of the large rotation rate on turbulence. Note that the 
s-K-e model gives the same results as for the no rotation case because it cannot 
account for the effect of rotation on the evolution of turbulence. 
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4. Backward-Facing Step Flows 
4.1 Numerical procedure 


For computational convenience, the non-dimensional form of the governing equa- 
tions is solved, in which 

X{ rr U* P 

< Xi > — — , — — , < P > rrj 1 

L ”t U "l pU "! (54) 

„ K £l„, V, ' ' 

< K >= Jjj-, < e >= — 3 ^-, < v t >= — — - — 

U; U? tf UrefL rt f 

where < > refers to a non-dimensional quantity, and L re j , and U Te f are the reference 
length and velocity, respectively. Accordingly, the flow Reynolds number is defined 
by 


Re = hllEl l l ( 55 ) 

v 

Hereafter, all the quantities will be of the non-dimensional form so that < > will be 
dropped for simplicity. 

In the steady-state and two dimensional cases (*1 = x,x 3 = y), the transport 
equations ( 1), ( 2), ( 48) and ( 49) can be written in the following general form 

[t^-(4- + + [v^ - (i + £•#,], = s * (56) 

He 0 $ tie 0 $ 

where <f> stands for 1 , U ( = U\)> V ( = U 2 ), K and c. For the momentum equations, 
the source term S# includes the cross-derivative diffusion and quadratic velocity 
gradient terms arising from equation ( 47). It can be seen that the non-dimensional 
equations are all of the same form as their dimensional counterparts, except that 
the kinematic molecular viscosity v is replaced by 1 / Re. 

The numerical method used to solve the system of equations ( 56) is a finite- 
volume procedure. It uses a non-staggered grid with all the dependent variables 
being stored at the geometric center of each control volume (Figure 3). The mo- 
mentum interpolation procedure of Rhie and Chow (1983) is used to avoid spurious 
oscillations usually associated with the non-staggered grid, and the pressure- velocity 
coupling is handled with the SIMPLEC algorithm (Van Doormal and Raithby, 1984). 
To ensure both accuracy and stability of numerical solution, the convection terms 
are approximated by a second-order accurate and bounded differencing scheme (Zhu, 
1991a), and all the other terms by the conventional central differencing scheme. As 
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a result, the discretized counterpart of equation ( 56) can be cast into the following 
linearized form 


E Ai4>c = Ai<fn + Sc (57) 

i 

where the coefficients Aj (l = W, E, S , N ), which relate the principal unknown <f>c 
to its neighbours <f>i (Figure 3), result from the discretization of the left-hand side 
terms of equation ( 56). The convection scheme used ensures that At > 0 so that 
the resulting coefficient matrix is always diagonally dominant. The strongly implicit 
procedure of Stone (1968) is used to solve the system of algebraic equations. The 
iterative solution process is considered converged when the maximum normalized 
residue of all the dependent variables is less than 10~ 4 . The details of the present 
numerical procedure are given in Rodi ef a 1. (1989) and Zhu (1991b). 

4.2 Numerical results 

The present model is then applied to the two backward-facing step flows exper- 
imentally studied by Kim, Kline and Johnston (1978) and Driver and Seegmiller 
(1985), from here on referred to as KKJ- and DS-cases, respectively. Figure 4 
shows the flow configuration and the Cartesian co-ordinate system used. Table 
1 gives the flow parameters for both cases; here the experimental reference free- 
stream velocity U re f and step height H, are taken as the reference quantities for 
non-dimensionalization. 


Table 1. Flow parameters 


case 

Re 

S 

L $ 

L e 

H. 

H d 

Uref 

DS 

37423 

1.5 

10 

40 

1 

8 

1 

KKJ 

44737 

0.6 

10 

40 

1 

2 

1 


Three types of boundaries are present, i.e. inlet, outlet and solid wall. At the 
inlet, the experimental data are available for the stream wise mean velocity U and 
the turbulent normal stresses uu and vv. K is calculated from these uu and vv with 
the assumption that 


and t by 


tutu = 


-(uu + vv) 


(58) 


C? / 4 KV 3 

e = -*-= , L = min(0.41Ay, 0.0855) (59) 
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where Ay is the distance from the wall and 8 is the boundary-layer thickness given 
in Table 1. At the outlet, the stream wise derivatives of the flow variables are set 
to zero. Influences of both inlet and outlet conditions on the solution are examined 
by changing the locations L, and L e , and it has been found that in both cases, 
the distances given in Table 1 are already sufficiently far away from the region of 
interest. In the earlier stage of this work, we tested several low Reynolds number 
K-e models including those of Chien (1982), Lam and Bremhorst (1981), Launder 
and Sharma (1974), Shih and Lumley (1992), and Yang and Shih (1992), but none 
of them was found to be able to yield satisfactory solutions. Similar findings were 
also reported in Avva et aJ. (1990), Shuen (1992) and So and Lai (1988). Therefore 
in this work, we use the standard wall function approach (Launder and Spalding, 
1974) to bridge the viscous sublayer near the wall. 

Two sets of non-uniform computational grids are used to examine the grid depen- 
dence of the solution; they contain 110x52 (coarse) and 199x91 (fine) points for the 
KKJ-case and 106x56 (coarse) and 201x109 (fine) points for the DS-case. Figures 
5(a) and 5(b) show the friction coefficient Cf at the bottom wall calculated with the 
s-K-c model and the present model; also included in figure 5(a) are the experimental 
data for the DS-case, but no such data are available for the KKJ-case. It can be 
seen that the grid refinement does produce some differences for the results of the 
present model, more noticeable in the KKJ-case, and this is also the case for the 
s-K-c results. This indicates that the solutions obtained on the coarse grids are not 
sufficiently close to the grid-independent stage. Recently, Thangam and Hur (1991) 
have conducted a highly-resolved calculation for the KKJ-case. They have found 
that quadrupling a 166x73 grid leads to only a minimal improvement. Therefore, 
the present results on the fine grids can be considered as grid-independent. For the 
DS-case, the fine grid computations with the s-K-e and present model required 703 
and 805 iterations, and took approximately 7.1 and 8.3 minutes of CPU time on the 
Cray YMP computer. In the following, only the fine grid results are presented. 

The wall friction coefficient Cf is a parameter that is very sensitive to the near- 
wall turbulence modeling. It is Cf that the various low Reynolds number K-c models 
tested predict much worse than those using wall functions. However, the influence 
of the near-wall turbulence modeling is only restricted in the near-wall regions. It 
is seen from figure 5(a) that both the s-K-c and present model largely underpredict 
the negative peak of Cf , pointing to limited accuracy of the wall function approach 
in the recirculation region. 

The computed and measured reattachment points are compared in Table 2. They 
are determined in the calculation from the point where C / goes to zero. Also included 
in Table 2 are the result of Obi et a 1. (1989) obtained with the Reynolds stress model 


17 



(RSM) and that of Sindir (1982) with a modified algebraic stress model (ASM). The 
reattachment point is a critical parameter which has often been used to assess the 
overall performance of turbulence models as well as numerical procedures. Table 2 
clearly demonstrates the significant improvement obtained with the present model. 
It is important to mention that this improvement is mainly due to the behavior of 
Cft in the present model, and that the anisotropic behavior of the turbulent stresses 
only makes a marginal contribution to it. 


Table 2. Comparison of reattachment points 


case 

experiment 

s-K-e present model 

RSM 

ASM 

DS 

6.1 

4.99 5.82 

- 

5.66 

KKJ 

7 ±0.5 

6.35 7.35 

6.44 

- 


Figures 6(a) and 6(b) show the comparison of computed and measured static 
pressure coefficient C p along the bottom wall. In both cases, the s-K-e is seen to 
predict premature pressure rises, which is consistent with its underprediction of the 
reattachment lengths, while the present model captures these pressure rises quite 
well. The good predictions of C p were reported in both works of Obi et a /. (1989) 
and Sindir (1982), using the RMS and ASM. The results of the present model are 
almost comparable to those of the RSM and ASM. Again, the improved predictions 
of C p are mainly attributed to the variation of C p . 

The streamwise mean velocity U profiles are shown in figures 7(a) and 7(b) at four 
different cross-sections. Here, the differences between the results of the s-K-e and 
present model are not substantial, as compared to other flow variables. The present 
model predicts reverse flows better than the s-K-e, but results in somewhat slower 
recovery in the vicinity of the reattachment point. Interestingly enough, such a slow 
recovery also exists in the RSM prediction of Obi et al. (1989). Further downstream, 
say at x=20 in figure 7(a), the results of the two models nearly coincide with each 
other. 

Finally, the comparisons of predicted and measured turbulent stresses u 2 , v 2 and 
St; are shown in figures 8 and 9 at various x-locations. In the KKJ-case, no ex- 
perimental data for the turbulent stresses are available in the recirculation region, 
and the reattachment point was found in the experiment to move forward and back- 
ward continuously around seven step heights downstream of the step, leaving an 
uncertainty of ±0.5 step height for the reattachment length. This also points to 
some uncertainty in the measured turbulent quantities in the recovery region. On 
the other hand, the experimental data in the DS-case should be considered more 
reliable because of the smaller uncertainty of the reattachment location, indicating 
a smaller unsteadiness of the flow. As compared with the s-K-e results in figures 
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8 and 9, it can be seen that the anisotropic terms increase u 2 while decreasing v 2 , 
leading to significant improvements in both u 7 and v 2 results. On the other hand, 
the anisotropic terms have little impact on the turbulent shear stress CtJ. These 
behaviours are clearly reflected in equations ( 41) which also hold qualitatively for 
the flows considered here. The improvement obtained by the present model in figure 
8 for uv is due to the reduction in (7 M . 


5. Conclusions 

A constitutive relation for the turbulent stresses has been derived by using invari- 
ance theory. The relation is valid only for turbulent flows at high Reynolds numbers 
because the influence of the molecular viscosity has not been taken into account in 
the analysis. Being a second rank tensor, the general form of the turbulent stress 
can be expressed as a series, in terms of the mean velocity gradients, of order up 
to 4, while the classical eddy- viscosity representation constitutes only a first-order 
approximation. For practical calculations, it may suffice to use the quadratic ap- 
proximation. The model coefficients are functions of the time scale ratios rj and 

which ensure the positiveness of the turbulent normal stresses - a realizability 
condition that most existing turbulence models are unable to satisfy. The present 
model has been applied to calculate the two different flows: rotating homogeneous 
shear flows and backward-facing step flows. The calculated results show that all flow 
variables are sensitive to the variation of C^, and that only the turbulent normal 
stresses are sensitive to the terms containing nonlinear mean velocity gradients. The 
values of the model coefficients given in this paper seem quite appropriate for both 
the test cases, but the value of C T i related to the cross-derivative quadratic terms 
has little impact on the flows considered. This indicates that C T \ may be further 
calibrated against other flows. The computed results have been compared in detail 
with the LES data for rotating homogeneous shear flows and the experimental data 
for backward-facing step flows. The comparisons show that the present model does 
provide significant improvement over the standard K-e model, and this improvement 
is achieved at an insignificant penalty to the computational efficiency and algorith- 
mic simplicity of the latter. The present model can also be expected to work well 
for simple inhomogeneous shear flows, as evidenced by its improved prediction in 
the region far downstream of the reattachment point where the flow tends to be of 
simple parabolic nature. 


19 



Acknowledgements 


The authors are grateful to Dr. Aamir Shabbir for his calculations of the rotat- 
ing homogeneous shear flows and for may helpful suggestins and discussions. The 
contribution of J.L. Lumley was supported in part by Contract No. AFOSR 89- 
0226, jointly funded by the U.S. Airforce Office of Scientific Research (Control and 
Aerospace Programs), and the U.S. Office of Naval Research, and in part by Grant 
No. F49620-92- J-0038, funded by the U.S. Airforce Office of Scientific Research 
(Aerospace Program). 


References 


1. R.K. Am, C.E. Smith and A.K. Singhal, 1990, “Comparative study of high 
and low Reynolds number versions of K-e models”, AIAA paper 90-0246. 

2. J. Bardina, J.H. Ferziger and W.C. Reynolds, 1983, “Improved turbulence 
models based on large-eddy simulation of homogeneous incompressible turbu- 
lent flows.” Rept. No.TF-19, Stanford University, Stanford, Ca. 

3. K.Y. Chien, 1982, “Predictions of channel and boundary-layer flows with a 
low- Reynolds- number turbulence model”, AIAA J., Vol.20, pp. 33-38. 

4. D.M. Driver and H.L. Seegmiller, 1985, “Features of a reattaching turbulent 
shear layer in divergent channel flow”, AIAA J ., Vol.23, pp. 163-171. 

5. P.G. Huang and M.A. Leschziner, 1985, “Stabilization of recirculating-flow 
computations performed with second moment closures and third order dis- 
cretization”, Proceedings of the 5th Symposium on Turbulent Shear Flows , 
Cornell University, pp. 5. 19-5.24. 

6. J. Kim, S.J. Kline and J.P. Johnston, 1978, “Investigation of separation and 
reattachment of a turbulent shear layer: Flow over a backward-facing step”, 
Rept. MD-37, Thermosciences Div., Dept, of Mech. Eng., Stanford Univer- 
sity. 

7. C.K.G. Lam and K. Bremhorst, 1981, “A modified form of K-e model for 
predicting wall turbulence”, J. Fluids Eng., Vol.103, pp.456-460. 


20 



8. B.E. Launder and B.I. Sharma, 1974, “Application of the energy-dissipation 
model of turbulence to the calculation of a flow near a spinning disk”, Letters 
in Heat and Mass transfer V ol.l, pp. 131-138. 

9. B.E. Launder and D.B. Spalding, 1974, “The numerical computation of tur- 
bulent flows”, Comput. Meths. App. Mech. Eng., Vol.3, pp. 269-289. 

10. M.A. Leschziner and W. Rodi, 1981, “Calculation of annular and twin parallel 
jets using various discretization schemes and turbulence model variations”, J. 
Fluids Eng., Vol.103, pp. 352-360. 

11. J.L. Lumley, 1970, “Toward a turbulent constitutive relation.” J. Fluid Mech., 
Vol.41, pp.413-434. 

12. J.L. Lumley, 1978, “Computational modeling of turbulent flows”, Adv. Appl. 
Mech., Vol.18, pp. 124-176. 

13. S. Obi, M. Peric and G. Scheuerer, 1989, “A finite-volume calculation pro- 
cedure for turbulent flows with second-order closure and co-located variable 
arrangement”, Report. LSTM 276/N/89, Lehrstuhl fur Stromungsmechanik, 
Universitat Erlangen- Nurnberg. 

14. R.S. Rivlin, 1955, “Further remarks on the stress deformation relations for 
isotropic materials”, J. Arch. Rati. Mech. Anal., Vol.4, pp.681-702. 

15. C.M. Rhie and W.L. Chow, 1983, “A numerical study of the turbulent flow past 
an isolated airfoil with trailing edge separation”, AIAA J., Vol.21, pp.1525- 
1532. 

16. W. Rodi, 1980, “Turbulence models and their application in hydraulics - A 
state of the art review”, Book Publication of the International Association for 
Hydraulic Research, Delft, the Netherlands. 

17. W. Rodi, S. Majumdar and B. Schonung, 1989, “Finite-volume method for two 
dimensional incompressible flows with complex boundaries”, Comput. Meths. 
App. Mech. Eng., Vol.75, pp. 369-392. 

18. R. Rubinstein and J.M. Barton, 1990, “Nonlinear Reynolds stress models and 
the renormalization group”, Phys. Fluids A 2, pp. 1472-1476. 

19. U. Schumann, 1977, “Realizability of Reynolds stress turbulence models”, 
Phys. Fluids, Vol.20, pp.721-725. 


21 



20. T.H. Shih, 1992, "Remarks on turbulent constitutive relations.” to appear in 
NASA TM. 

21. T.H. Shih and J.L. Lumley, 1992, “Kolmogorov behavior of near-wall turbu- 
lence and its application in turbulence modeling”. NASA TM 105663, also in 
Int. J. Comput. Fluid Dynamics, Vol.l. 

22. J.S. Shuen, 1992, Private communication. 

23. M. Sindir, 1982, “Numerical study of separating and reattaching flows in a 
backward-facing step geometry”, Ph.D. Thesis, University of California at 
Davis. 

24. R.M.C. So and Y.G. Lai, 1988, “Low- Reynolds-number modelling of flows over 
a backward-facing step”, J. Appl. Math. Phys. (ZAMP), Vol.39, pp. 13-27. 

25. C.G. Speziale, 1987, “On nonlinear K-l and K-e models of turbulence”, J. 
Fluid Mech., Vol.178, pp.459-475. 

26. C.G. Speziale and S. Thangam, 1992, “Analysis of an RNG based turbulence 
model for separated flows”, NASA CR- 189600, ICASE Report. No.92-3. 

27. H.L. Stone, 1968, “Iterative solution of implicit approximations of multidimen- 
sional partial differential equation”, SIAM J. Num. Anal., Vol.5, pp. 530-558. 

28. D.B. Taulbee, 1992, “An improved algebraic Reynolds stress model and corre- 
sponding nonlinear stress model”, Phys. Fluids A, Vol.4, No. 11, pp.2555-2561. 

29. S. Thangam and N. Hur, 1991, “A highly-resolved numerical study of turbulent 
separated flow past a backward-facing step”, Int. J. Eng. Sci., Vol.29, pp.607- 
615. 

30. J.P. Van Doormal and G.D. Raithby, 1984, “Enhancements of the SIMPLE 
method for predicting incompressible fluid flows”, Num. Heat Trans., Vol.7, 
pp.147-163. 

31. Z. Yang and T.H. Shih, 1992, “A new time scale based K-e model for near 
wall turbulence”, NASA TM 105768. 

32. A. Yoshizawa, 1984, “Statistical analysis of the deriation of the Reynolds stress 
from its eddy-viscosity representation”, Phys. Fluids, Vol.27, pp. 1377-1387. 

33. J. Zhu, 1991a, “A low diffusive and oscillation-free convection scheme”, 
Comm, App. Num. Methods., Vol.7, pp.225-232. 


22 



34. J. Zhu, 1991b, “FAST-2D: A computer program for numerical simulation of 
two-dimensional incompressible flows with complex boundaries”, 

Rept. No. 690, Institute for Hydromechanics, University of Karlsruhe. 


23 



Appendix A 


Generalized Cayley-Hamilton Formulas 

RivEn (1955) showed that there are several generalized Cayley-Hamilton formulas 
relating matrices (product of several matrices A,B t C ...) of higher extension to 
matrices of lower extension. Some of them are Ested here for latter use. 

ABC + ACB + BCA + BAC + CAB + CBA - A(trBC - trB trC) 

— B(trCA — trC trA ) — C(trAB — trA trB) 

-{BC + CA)trA - ( CA + AC)trB - ( AB + BA)trC (A.l) 

-I(trA trB trC - trA trBC - trB trCA 
—trC trAB + trABC + trCBA) = 0 
Replacing C with A and B in Eq.( A.l), respectively, we obtain 
ABA = -A 2 B - BA 2 + A{trAB - trA trB) 

+ i B(trA 2 - trA trA) + (AB + BA)trA + A 2 trB ( A 2 ) 

-\-I[trA 2 B — trA trAB 4 - ^ trB(trA trA — trA 2 )\ 

and 

BAB = -B 2 A - AB 2 + B(trBA - trB trA) 

-f ~A(trB 2 - trB trA) + (BA + AB)trB + B 2 trA ( A 3 ) 

+I[trB 2 A - trB trBA-\-\trA(trB trB - trB 2 )) 

which indicate that the matrices ABA and BAB of extension 3 can be expressed 
by polynomials of matrices of extension 2 or less. 
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Multiplying Eq.( A. 2) from the left and the right by A and Eq.( A. 3) by B , and 
adding them correspondingly, we obtain following two relations: 

ABA 2 + A 2 BA = ABA trA + A 2 trAB + A(trA 2 B - trA trAB) 

(A.4) 

— B detA + I detA trB 


and 


BAB 2 + B 2 AB = BAB trB + B 2 trB A + B(trB 2 A - trB trB A) 
—A detB + I detB trA 


(A.5) 


Replacing B with B 2 in Eq.( A.4) and A with A 2 in Eq.( A.5) give 
AB 2 A 2 + A 2 B 2 A = AB 2 A trA + A 2 trAB 2 
+A(trA 2 B 2 - trA trAB 2 ) - B 2 detA + 1 detA trB 2 


(A.6) 


BA 2 B 2 + B 2 A 2 B = BA 2 B trB + B 2 trB A 2 

+B(trB 2 A 2 - trB trB A 2 ) - A 2 detB + I detB trA 2 
Replacing B with B 2 in Eq.( A. 2) and A with A 2 in Eq.( A.3) yield 

AB 2 A = - A 2 B 2 - B 2 A 2 + A(trAB 2 - trA trB 2 ) 

+ i B 2 {trA 2 - trA trA) + (AB 2 + B 2 A)trA + A 2 trB 2 
2 

+I[trA 2 B 2 - trA trAB 2 + \trB 2 {trA trA - trA 2 )] 

and 

BA 2 B = -B 2 A 2 - A 2 B 2 + B(trBA 2 - trB trA 2 ) 

+- A 2 (trB 2 - trB trA 2 ) + (BA 2 + A 2 B)trB + B 2 trA 2 
2 

+I[trB 2 A 2 - trB trB A 2 + \ trA 2 (trB trB - trB 2 )] 


(A.T) 


(A.8) 


(A.9) 


Eqs.( A.8) and ( A.9) indicate that the matrices AB 2 A and BA 2 B of extension 3 can 
be expressed by polynomials of matrices of extension 2 or less. Therefore, the right 
hand sides of Eqs.( A.6) and ( A. 7) are also polynomials of matrices of extension of 
2 or less. 
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Appendix B 


Number of Independent Tensors Formed by Two Tensors 


Let us show that the number of independent tensors formed with two general 
tensor A and B is 18. 

Rivlin (1955) showed that any matrix product in two 3x3 matrices may be 
expressed as a polynomial in these matrices of extension 4 or less. Suppose we have 
a matrix product II of extension 5: 


H = ABA 2 B 2 A 


(B.l) 


This can be written as 

n = AC A (B.2) 

where C = BA 2 B 2 . From Eq.( A.2), II may be viewed as a polynomial of matrices 
in A and C of extension 2 or less. C itself is a matrix in A and B of extension 
3 so that II may be expressed by a polynomial in A and B of extension 4 or less. 
Therefore, we only need to consider the possible tensors of extension 4 or less formed 
by A and B. 

We may show that there are only two independent tensors of extension 4. The 
possible tensors of extension 4 are the following 8 tensors: 


ABA 2 B 2 , BAB 2 A 2 , A 2 BAB 2 y B 2 ABA 
AB 2 A 2 B , BA 2 B 2 A , A 2 B 2 AB , B 2 A 2 BA. 


(B.3) 


With Eq.( A.4), A 2 BAB 2 can be expressed by ABA 2 B 2 + • • Similarly, with 
Eq.( A.5), B 2 ABA 2 = -BAB 2 A 2 + with Eq.( A.7), AB 2 A 2 B = -ABA 2 B 2 

+ •••; with Eq.( A.6), BA 2 B 2 A = -BAB 2 A 2 + •••; with Eqs.( A.5) and ( A.4), 
A?B 2 AB = ABA 2 B 2 + • • •; with Eqs.( A.4) and ( A.5), B 2 A 2 BA = BAB 2 A 2 + • • •; 
where • • • represents a polynomial in A and B of extension 3 or less. As a result, 
only two tensors of extension 4 in Eq.( B.3) are independent, and we select them as 


ABA 2 B 2 , BAB 2 A 2 


(BA) 


26 


Ubt-WAL PAGE lb 

Of POOR QUALITY 



Now we show that there are only four independent tensors of extension 3. The 
possible tensors of extension 3 are the following 8 tensors: 

ABA 3 , A 2 BA , BAB 3 , B 2 AB, AB 2 A 3 , A 3 B 3 A, BA 3 B \ B 3 A 3 B. (B.5) 

Using Eqs.( A.4), ( A. 5), ( A.6) and ( A. 7), we find that only four of them are 
independent. Let us select them as 

ABA 2 , BAB 2 , AB 3 A 3 , BA 3 B 3 . (B.6) 

Furthermore, there are eight independent tensors of extension 2: 

AB, BA, AB 3 , B 2 A, A 3 B, BA 3 , A 3 B 3 , B 3 A 3 (B.7) 

and four independent tensors of extension 1: 

A, A 2 , B, B 2 . (B.8) 

Therefore, we have proved that only 18 tensors can be formed independently by 
two general tensors. 
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Figure 2. Evolution of turbulent kinetic energy with time. 
: s-K-e; : present model; •: LES. 
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Figure 3. Typical control volume centered at C and related nodes 



Figure 4. Backward-facing step geometry 
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Figure 5. Friction coefficient C f along the bottom wall. 

: s-K-e, fine grid; : present model, fine grid; 

: present model, coarse grid; •: experiment. 
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Figure 7. Stream wise mean velocity {/-profiles. 

• s-K-e; : present model; •: experiment 
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Figure 9. Turbulent stress profiles in KKJ-case. 

- • s-K-e; : present model; •: experiment 
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